Introduction

      Benthic macroinvertebrate communities show structural and functional differences depending on natural landscape or stream characteristics such as stream size, geology, climate, and elevation (Vannote et al. 1980; Dodds et al., 2015; Cushing et al., 2006). These landscape differences determine the water chemistry, food resources and habitat available for macroinvertebrates. Certain taxa are uniquely adapted to withstand conditions that occur across landscapes naturally, such as organisms that are adapted to cold water or have physical properties to protect them from abrasion in sandier habitats (Poff et al., 2006). In the absence of extreme human modifications, geology, topography, climate, and stream size result in somewhat predictable macroinvertebrate communities that reflect natural disturbance regimes (Dodds et al., 2015; Poff et al., 2006). While considering variation across benthic communities that could be attributed to stress, it is important to understand and recognize the biotic differences that occur naturally.

      One method to identify which natural landscape and stream characteristics are having a measurable effect on the organisms found in a macroinvertebrate community is to use multivariate statistical techniques such as Non-metric Multidimensional Scaling (NMS) (McCune and Grace, 2002). Reference taxa data can be graphically reviewed using Non-metric Multidimensional Scaling (NMS) ordination and statistically explored for differences using the Multiple Response Permutation Procedure (MRPP) (McCune and Grace, 2002). By overlaying the stations natural landscape or stream characteristics, the patterns that may be associated with differences in natural characteristics emerge.

      Identifying the important natural landscape and stream characteristics to a biotic community is an essential first step for metric selection and Index of Biotic Integrity (IBI) development. If benthic communities are variable due to natural characteristics, then associating biological communities with specific stressors will be inaccurate. The Virginia Department of Environmental Quality is undertaking the development of a genus-level IBI to replace the family level Virginia Stream Condition Index (TetraTech, 2003; VDEQ, 2006). Virginia has a wide variety of topographic and geologic characteristics and includes both mountainous and coastal streams. Within Virginia there are seven Level 3 ecoregions and three coarser physiographic regions (called bioregions herein). Three IBI’s were developed and are currently being used for the family level metrics, including a non-coastal metric (VSCI), coastal for ecoregion 65 (Virginia Coastal Plain Macroinvertebrate Index (VCPMI) - Chowan), and coastal for ecoregion 63 (VCPMI + Chowan)(VDEQ, 2013). As a first step in genus metric development, VADEQ performed a series of NMS for reference benthic communities to identify clusters that could be attributed to natural landscape characteristics and may require individual metric development.

Methods

Study Design

      Reference stations were identified using a series of landscape and water quality filters including VSCI scores that were above the impairment threshold. Across Virginia and part of West Virginia (Central Applications Ecoregion), 1045 stations were identified as reference stations after evaluation from the regional biologists (Ref stress Map-https://rconnect.deq.virginia.gov/RefStressMap/). Some of West Virginia reference data was used to increase the number of stations located in the Central Appalachians ecoregion. The macroinvertebrate collection, sorting, and identification methods used are similar to Virginia; therefore, it was deemed that this data was appropriate for usage. The dataset resulted in 696 stressed samples (n Sites= 233), 1129 intermediately stressed samples (n sites=417), 858 Reference samples (n sites=336 ), 100 Reference Piedmont samples (n sites=40), and 117 Reference West Virginia samples (n sites= 61). Reference, Reference Piedmont and West Virginia samples were all used in the analysis below.

      A macroinvertebrate community abundance matrix was created from the reference stations identified. All taxa denoted as excluded by the biologists during identification were not used in this evaluation. Taxa that were identified to Order were excluded from the matrix to reduce noise in the dataset. If stations appeared in the dataset more than 4 times, then the most recent four samples were used and the remaining data was removed. Any samples that had less than 100 total taxa in a sample were also removed. Taxa that occurred in the dataset less than or equal to 5% of the time were removed from the dataset to eliminate rare taxa. All abundance data was then log_10+1 transformed for use in the ordinations.

      A small review group was formed of EPA and VADEQ experts in IBI development to review and provide input on the metric development process. Based on feedback from the group, iterations of the NMS technique were conducted to evaluate the removal of 1%, 3%, and 10% of rare taxa and the removal of Chironomidae, which were ubiquitous in all samples. The suggested changes did not appear to significantly change the outcome of the NMS ordinations; therefore they are not included in the results section below.

Statistical Analysis

      Non-metric Multiple Dimensional Scaling (NMS) was used to examine the reference dataset for patterns associated with landscape or stream characteristics. This ordination technique plots macroinvertebrate community abundance data in 3-D space to identify if certain station’s communities cluster together in a non-random way. Distances among samples were calculated using a Bray-Curtis Dissimilarity distance. Rather than plotting on just two axes (x and y), NMS allows the user to specify the number of dimensions, allowing the user to see more refined clustering graphically. We evaluated 3 axes; however, axes 1 and 2 accounted for the most variation and are shown throughout this document. NMS will iteratively seek the best solution for the dataset and stop once it has found the best solution or at a pre-determined number of iterations. The best solution is considered when the stress value is at its lowest. To evaluate the strength of the NMS stress, we used Clarke’s method that states if the stress is at or above 0.20 then the fit is not different than a random ordination. Stress below 0.20 is considered different than random and provides an ordination where inferences can be drawn (McCune and Grace 2002- Chapter 16).

      Multiple Response Permutation Procedure (MRPP) provides a test to determine if there is a significant difference between two or more groups. This is a non-parametric test and is similar to tests like the Analysis of Variance (ANOVA) (Biodini et al., 1988). Bray Curtis distance was used for this statistical test. The test statistic, A describes the within-group similarity compared to a random expectation. Within community ecology data, a test statistic of greater than 0.3 is considered fairly high and indicates that groups are considered different from one another (McCune & Grace, 2002). An A statistic of 1 indicates that all within group communities are identical. The NMS ordinations and MRPP were conducted in R version 4.0.3 (2020-10-10) – “Bunny-Wunnies Freak Out” using the vegan package (version 2.5-7).

Results

All Reference data

      The results of the initial NMS ordinations using all reference data indicated that several landscape and stream characteristics did not demonstrate a difference in benthic macroinvertebrate communities. These characteristics include season, sample method, and trout waters (Table 1, A<0.07).

      The NMDS and MRPP results did, however, indicated that there are differences in benthic reference communities based on Bioregion and Ecoregion. The Coastal bioregion visibly separated from the Piedmont and Mountain bioregions and the Mid-Atlantic Coastal Plain and Southeastern Piedmont ecoregions were distinct from the other ecoregions (Figure 1 and 2).This relationship was evidenced by the MRRP analysis that showed relatively high A-statistics for these parameters (Bioregion= 0.0714 and Ecoregion=0.0970).Other A -statistics that were relatively high were likely correlated with bioregion and do not represent a meaningful difference in biological communities (Admin Region, Sample Method, Water Quality Standards). Because of the observed differences in benthic communities in Coastal and non-coastal systems the decision was made to separate the coastal and non-coastal systems for further ordination analyses. This decision is supported by biologists’ observations and expert group panel discussions and is consistent with the development of Virginia’s family level IBI (VDEQ, 2013).

Figure 1: Bioregion displayed on NMDS with all reference sites

Figure 1: Bioregion displayed on NMDS with all reference sites

Figure 2: Ecoregion displayed for all reference sites.

Figure 2: Ecoregion displayed for all reference sites.

Coastal Reference Streams NMDS

      The coastal reference streams (Ecoregions Mid-Atlantic Coastal Plain and Southeastern Plains) were separated from the non-coastal streams and reevaluated to identify any patterns that exist within coastal streams. There was minimal separation between spring and fall samples. Ecoregion separation existed between the Middle Atlantic Coast and Southeastern Plains. Some of this separation can be attributed to different basins (Figure 5). However, the coastal reference dataset is limited and samples may be geographically clustered where there is not adequate coverage to create separate IBIs for ecoregions (Figure 3 and 4). MRPP results did not support further separation of groups (Table 2).

Figure 3: Season displayed on NMS with coastal data only.

Figure 3: Season displayed on NMS with coastal data only.

Figure 4: Ecoregion displayed on NMS with coastal data only.

Figure 4: Ecoregion displayed on NMS with coastal data only.

Figure 5: Basin displayed on NMS with coastal data only.

Figure 5: Basin displayed on NMS with coastal data only.

Non-coastal Reference Streams NMDS

      Reference streams located in the Mid-Atlantic Coastal Plain ecoregion and Southeastern Piedmont ecoregion were removed from the analysis below. Rare taxa were removed (<=5%) from the community data matrix. After the removal of Coastal streams, several of the environmental variables evaluated did not show differences in macroinvertebrate communities including Administrative office, Basin, Sample Method, Water quality standards, natural trout waters, and ecoregion (Figure 6-10 & Table 3).

Figure 6: Administrative region displayed on non-coastal reference stream ordination.

Figure 6: Administrative region displayed on non-coastal reference stream ordination.

Figure 7: Stream basin displayed on non-coastal reference stream ordinations.

Figure 7: Stream basin displayed on non-coastal reference stream ordinations.

Figure 8: Sample Method displayed on non-coastal reference stream ordinations.

Figure 8: Sample Method displayed on non-coastal reference stream ordinations.

Figure 9: Non-coastal reference stream ordinations with water quality standard groupings.

Figure 9: Non-coastal reference stream ordinations with water quality standard groupings.

Figure 10: Non-coastal reference stream ordinations with natural trout water designation groupings.

Figure 10: Non-coastal reference stream ordinations with natural trout water designation groupings.

      Patterns were visible between several environmental parameters when coastal streams were removed and warrant additional investigation during the development of a Virginia IBI. These include Season, ecoregion, bioregion, and stream order (Figure 11-15). Stream order was categorized as small (Strahler Order 1-3) or large (Strahler Order 4-6) for better visualization. The ellipses shown on the figures below represent 90% of the data by category.

Figure 11: Non-coastal reference stream ordinations with Seasonal groupings.

Figure 11: Non-coastal reference stream ordinations with Seasonal groupings.

Figure 12: Non-coastal reference stream ordinations with bioregion groupings.

Figure 12: Non-coastal reference stream ordinations with bioregion groupings.

Figure 13: Non-coastal reference stream ordinations with Stream Order groupings.

Figure 13: Non-coastal reference stream ordinations with Stream Order groupings.

Figure 14: Non-coastal reference stream ordinations with ecoregion groupings.

Figure 14: Non-coastal reference stream ordinations with ecoregion groupings.

Figure 15: Non-coastal reference stream ordination with bioregion, stream order, and season category groupings.

Figure 15: Non-coastal reference stream ordination with bioregion, stream order, and season category groupings.

Season Separation

      Based on the visual differences between fall and spring seasons observed in among the non-coastal stations, it was deemed necessary to separate the seasons to better evaluate patterns in stream order, ecoregion, and bioregion (Figure 16-22)(Table 6). This evaluation will help us determine if further separation by environmental characteristics is necessary for Virginia IBI development.

Figure 16: Non-coastastal spring reference streams ordinations with ecoregion groupings

Figure 16: Non-coastastal spring reference streams ordinations with ecoregion groupings

Figure 17: Non-coastal fall reference stream ordinations with ecoregion groupings.

Figure 17: Non-coastal fall reference stream ordinations with ecoregion groupings.

Figure 18: Non-coastal spring reference stream ordinations with bioregion groupings.

Figure 18: Non-coastal spring reference stream ordinations with bioregion groupings.

Figure 19: Non-coastal fall reference stream ordination with bioregion groupings.

Figure 19: Non-coastal fall reference stream ordination with bioregion groupings.

Figure 20: Non-coastal spring reference stream ordination with stream order cateogry groupings.

Figure 20: Non-coastal spring reference stream ordination with stream order cateogry groupings.

Figure 21: Non-coastal fall reference stream ordination with stream order category groupings.

Figure 21: Non-coastal fall reference stream ordination with stream order category groupings.

Recommendations

      Based on the NMS ordinations shown in the results section above, we recommend several environmental variables be investigated further for differences in benthic macroinvertebrate community structure during IBI development. This investigative step is important in metric development so that the metric accurately represents the quality of a stream and minimizes variation caused by environmental characteristics. Environmental variation in benthic macroinvertebrate communities is inevitable; however, the ability to identify and minimize the environmental noise will allow biologists to better identify and understand stressed conditions. We recommend further investigation and separation of the following environmental parameters during the next phase of genus metric development (n represents the number of reference sites in each category):

Conclusion

      Several environmental characteristics in Virginia tend to have a significant effect on the macroinvertebrate community found in reference streams. These characteristics must be accounted for in order to create the most accurate genus level IBI. Specifically, the ordinations show clear separation in reference communities found in coastal systems and non-coastal systems. This separation makes sense given the macroinvertebrate sampling method (riffle vs. MACS) and the environmental variables often associated with coastal streams versus non-coastal streams, including slower moving flow, finer substrate, and potentially brackish water. These environmental variables directly affect the habitat and water quality resulting in lower dissolved oxygen concentrations, higher conductivity levels, and higher pH. The resulting macroinvertebrate community will have adaptations to withstand these conditions.

      Among coastal streams, no other environmental characteristics were identified as showing separation in the benthic macroinvertebrate community. However, among the non-coastal streams there was slight separation in the ordinations between spring and fall communities. Spring sampling is conducted between March and June and fall sampling is conducted between September and November. In Virginia, spring and fall seasons can vary greatly but Spring is often characterized by higher precipitation and cooler temperatures while Fall is often drier and warmer. There was also separation between ecoregions in coastal streams. Further investigation is warranted to identify the need for two separate IBI to characterize spring and fall macroinvertebrate communities.

      Non-coastal spring and fall communities were separated to analyze patterns in other environmental variables among season. Bioregion and stream order showed some minor grouping between categories for both spring and fall. Further investigation is warranted to identify if these categories require their own IBI. However, the number of streams in each category is reduced due to the number of categories required and may limit the ability to develop a robust IBI for each group.

References

Biondini, Mario E., Paul W. Mielke, and Kenneth J. Berry. “Data-dependent permutation techniques for the analysis of ecological data.” Vegetatio 75.3 (1988): 161-168.

Dodds, W. K., Gido, K., Whiles, M. R., Daniels, M. D., & Grudzinski, B. P. (2015). The stream biome gradient concept: Factors controlling lotic systems across broad biogeographic scales. Freshwater Science, 34(1), 1-19.

McCune, B., & Grace, J. B. (2002.) Analysis of ecological communities. Gleneden Beach, Oregon: MjM Software Design. Chapter 24.

Poff, N. L., Olden, J. D., Vieira, N. K., Finn, D. S., Simmons, M. P., & Kondratieff, B. C. (2006). Functional trait niches of North American lotic insects: traits-based ecological applications in light of phylogenetic relationships. Journal of the North American Benthological Society, 25(4), 730-755.

R Core Team (2019). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

Tetra Tech, Inc. 2003. A Stream Condition Index for Virginia Non-Coastal Streams. Tetra Tech, Inc. Owings Mills, Maryland. Prepared for Virginia Department of Environmental Quality, Richmond, Virginia. https://www.deq.virginia.gov/home/showpublisheddocument/4317/637461491373170000

Vannote, R. L., Minshall, G. W., Cummins, K. W., Sedell, J. R., & Cushing, C. E. (1980). The river continuum concept. Canadian journal of fisheries and aquatic sciences, 37(1), 130-137.

Virginia Department of Environmental Quality. 2006. Using Probabilistic Monitoring Data to Validate the Non-Coastal Virginia Stream Condition Index. Technical Bulletin. WQA/2006-001.Richmond, Virginia. https://www.deq.virginia.gov/home/showpublisheddocument/4319/637461491379900000

Virginia Department of Environmental Quality. 2013. The Virginia Coastal Plain Macroinvertebrate Index. Technical Bulletin. WQA/2013-002.Richmond, Virginia. https://www.deq.virginia.gov/home/showpublisheddocument/4315/637461491365370000